home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-04 / bipl.zip / PROCS.ZIP / MATH.ICN < prev    next >
Text File  |  1992-09-28  |  7KB  |  239 lines
  1. ############################################################################
  2. #
  3. #    File:     math.icn
  4. #
  5. #    Subject:  Procedures to perform mathematical computations
  6. #
  7. #    Author:   George D. Yee
  8. #
  9. #    Date:     June 10, 1988
  10. #
  11. ###########################################################################
  12. #
  13. #  Note:
  14. #     Version 8 of Icon supports most of the procedures that follow as
  15. #  built-in functions. The procedures here should not be used unless the
  16. #  corresponding functions are disabled.
  17. #
  18. ############################################################################
  19. #  
  20. #  The following procedures compute standard trigonometric func-
  21. #  tions.  The arguments are in radians.
  22. #  
  23. #       sin(x)      sine of x
  24. #  
  25. #       cos(x)      cosine of x
  26. #  
  27. #       tan(x)      tangent of x
  28. #  
  29. #       asin(x)     arc sine of x in the range -pi/2 to pi/2
  30. #  
  31. #       acos(x)     arc cosine of x in the range 0 to pi
  32. #  
  33. #       atan(x)     arc tangent of x in the range -pi/2 to pi/2
  34. #  
  35. #       atan2(y,x)  arc tangent of x/y in the range -pi to pi
  36. #  
  37. #  The following procedures convert from degrees to radians and con-
  38. #  versely:
  39. #  
  40. #       dtor(d)     radian equivalent of d
  41. #  
  42. #       rtod(r)     degree equivalent of r
  43. #  
  44. #  The following additional procedures are available:
  45. #  
  46. #       sqrt(x)     square root of x
  47. #  
  48. #       exp(x)      exponential function of x
  49. #  
  50. #       log(x)      natural logarithm of x
  51. #  
  52. #       log10(x)    base-10 logarithm of x
  53. #  
  54. #       floor(x)    largest integer not greater than x
  55. #  
  56. #       ceil(x)     smallest integer nor less than x
  57. #  
  58. #  Failure Conditions: asin(x) and acos(x) fail if the absolute
  59. #  value of x is greater than one. sqrt(x), log(x), and log10(x)
  60. #  fail if x is less than zero.
  61. #  
  62. ############################################################################
  63.  
  64. procedure sin(x)
  65.    return _sinus(numeric(x),0)
  66. end
  67.  
  68. procedure cos(x)
  69.    return _sinus(abs(numeric(x)),1)
  70. end
  71.  
  72. procedure tan(x)
  73.    return sin(x) / (0.0 ~= cos(x))
  74. end
  75.  
  76. # atan returns the value of the arctangent of its
  77. # argument in the range [-pi/2,pi/2].
  78. procedure atan(x)
  79.    if numeric(x) then
  80.       return if x > 0.0 then _satan(x) else -_satan(-x)
  81. end
  82.  
  83. # atan2 returns the arctangent of y/x
  84. # in the range [-pi,pi].
  85. procedure atan2(y,x)
  86.    local r
  87.    static pi
  88.    initial pi := 3.141592653589793238462643
  89.    return if numeric(y) & numeric(x) then {
  90.       if x > 0.0 then
  91.          atan(y/x)
  92.       else if x < 0.0 then {
  93.          r := pi - atan(abs(y/x))
  94.          if y >= 0.0 then r else -r
  95.          }
  96.       else if x = y = 0.0 then
  97.          0.0         # special value if both x and y are zero
  98.       else
  99.          if y >= 0.0 then pi/2.0 else -pi/2.0
  100.       }
  101. end
  102.  
  103. procedure asin(x)
  104.    if abs(numeric(x)) <= 1.0 then
  105.       return atan2(x, (1.0-(x^2))^0.5)
  106. end
  107.  
  108. procedure acos(x)
  109.    return 1.570796326794896619231e0 - asin(x)
  110. end
  111.  
  112. procedure dtor(deg)
  113.    return numeric(deg)/57.29577951308232
  114. end
  115.  
  116. procedure rtod(rad)
  117.    return numeric(rad)*57.29577951308232
  118. end
  119.  
  120. procedure sqrt(x)
  121.     return (0.0 <= numeric(x)) ^ 0.5
  122. end
  123.  
  124. procedure floor(x)
  125.    return if numeric(x) then
  126.       if x>=0.0 | real(x)=integer(x) then integer(x) else -integer(-x+1)
  127. end
  128.  
  129. procedure ceil(x)
  130.    return -floor(-numeric(x))
  131. end
  132.  
  133. procedure log(x)
  134.    local z, zsq, ex
  135.    static log2, sqrto2, p0, p1, p2, p3, q0, q1, q2
  136.    initial {
  137.       # The coefficients are #2705 from Hart & Cheney. (19.38D)
  138.       log2   :=  0.693147180559945309e0
  139.       sqrto2 :=  0.707106781186547524e0
  140.       p0     := -0.240139179559210510e2
  141.       p1     :=  0.309572928215376501e2
  142.       p2     := -0.963769093368686593e1
  143.       p3     :=  0.421087371217979714e0
  144.       q0     := -0.120069589779605255e2
  145.       q1     :=  0.194809660700889731e2
  146.       q2     := -0.891110902798312337e1
  147.       }
  148.    if numeric(x) > 0.0 then {
  149.       ex := 0
  150.       while x >= 1.0 do {
  151.          x /:= 2.0
  152.          ex +:= 1
  153.          }
  154.       while x < 0.5 do {
  155.          x *:= 2.0
  156.          ex -:= 1
  157.          }
  158.       if x < sqrto2 then {
  159.          x *:= 2.0
  160.          ex -:= 1
  161.          }
  162.       return ((((p3*(zsq:=(z:=(x-1.0)/(x+1.0))^2)+p2)*zsq+p1)*zsq+p0)/
  163.              (((1.0*zsq+q2)*zsq+q1)*zsq+q0))*z+ex*log2
  164.       }
  165. end
  166.  
  167. procedure exp(x)
  168.    return 2.718281828459045235360287 ^ numeric(x)
  169. end
  170.  
  171. procedure log10(x)
  172.    return log(x)/2.30258509299404568402
  173. end
  174.  
  175. procedure _sinus(x,quad)
  176.    local ysq, y, k
  177.    static twoopi, p0, p1, p2, p3, p4, q0, q1, q2, q3
  178.    initial {
  179.       # Coefficients are #3370 from Hart & Cheney (18.80D).
  180.       twoopi :=  0.63661977236758134308
  181.       p0     :=  0.1357884097877375669092680e8
  182.       p1     := -0.4942908100902844161158627e7
  183.       p2     :=  0.4401030535375266501944918e6
  184.       p3     := -0.1384727249982452873054457e5
  185.       p4     :=  0.1459688406665768722226959e3
  186.       q0     :=  0.8644558652922534429915149e7
  187.       q1     :=  0.4081792252343299749395779e6
  188.       q2     :=  0.9463096101538208180571257e4
  189.       q3     :=  0.1326534908786136358911494e3
  190.       }
  191.    if x < 0.0 then {
  192.       x := -x
  193.       quad +:= 2
  194.       }
  195.    y := (x *:= twoopi) - (k := integer(x))
  196.    if (quad := (quad + k) % 4) = (1|3) then
  197.       y := 1.0 - y
  198.    if quad > 1 then
  199.       y := -y
  200.    return (((((p4*(ysq:=y^2)+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y) /
  201.            ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0)
  202. end
  203.  
  204. procedure _satan(x)
  205.    static sq2p1,sq2m1,pio2,pio4
  206.    initial {
  207.       sq2p1 := 2.414213562373095048802e0
  208.       sq2m1 := 0.414213562373095048802e0
  209.       pio2  := 1.570796326794896619231e0
  210.       pio4  := 0.785398163397448309615e0
  211.       }
  212.    return if x < sq2m1 then
  213.              _xatan(x)
  214.           else if x > sq2p1 then
  215.              pio2 - _xatan(1.0/x)
  216.           else
  217.              pio4 + _xatan((x-1.0)/(x+1.0))
  218. end
  219.  
  220. procedure _xatan(x)
  221.    local xsq
  222.    static p4,p3,p2,p1,p0,q4,q3,q2,q1,q0
  223.    initial {
  224.       # coefficients are #5077 from Hart & Cheney. (19.56D)
  225.       p4    := 0.161536412982230228262e2
  226.       p3    := 0.26842548195503973794141e3
  227.       p2    := 0.11530293515404850115428136e4
  228.       p1    := 0.178040631643319697105464587e4
  229.       p0    := 0.89678597403663861959987488e3
  230.       q4    := 0.5895697050844462222791e2
  231.       q3    := 0.536265374031215315104235e3
  232.       q2    := 0.16667838148816337184521798e4
  233.       q1    := 0.207933497444540981287275926e4
  234.       q0    := 0.89678597403663861962481162e3
  235.       }
  236.    return x * ((((p4*(xsq:=x^2)+p3)*xsq+p2)*xsq+p1)*xsq+p0) /
  237.           (((((xsq+q4)*xsq+q3)*xsq+q2)*xsq+q1)*xsq+q0)
  238. end
  239.